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ABSTRACT 


A  large  sample  theory  is  presented  for  birth  and  death  queueing  processes 
which  are  ergodic  and  metrically  transitive.  The  theory  is  applied  to 
make  inferences  about  how  arrival  and  service  rates  vary  with  the  number 
in  the  system.  Likelihood  ratio  tests  and  maximum  likelihood  estimators 
are  derived  for  simple  models  which  describe  this  variation.  Composite 
hypotheses  such  as  that  the  arrival  rate  does  not  vary  with  the  number 
in  the  system  are  considered.  A  numerical  example  illustrating  these 
results  is  presented  which  includes  the  testing  of  both  true  and  false 
composite  hypotheses. 


PROBLEMS  OF  STATISTICAL  INFERENCE  FOR  BIRTH  AND  DEATH 
QUEUEING  MODELS 


I.  Introduction 

The  orientation  in  queueing  theory  is  to  start  with  a  model  and  determine 
the  behavior  of  a  system  from  it.  In  applications,  however,  one  has  the  problem 
of  deciding  what  model  is  appropriate  from  observing  the  behavior  of  a  system. 

For  queueing  processes,  the  problem  has  been  untouched  in  the  literature 
except  for  Clark  [3],  while  for  Markov  processes,  much  has  been  done.  A  large 
sample  theory  based  on  mwviimim  likelihood  for  Markov  processes  has  been  developed 
in  a  monograph  by  Billingsley  [l],  while  his  paper  [2]  contains  a  large  biblio¬ 
graphy.  The  results  parallel  the  classical  large  sample  theory  as  presented  in 
Wilks  [6],  Birth  and  death  queueing  processes  are  also  Markov  processes,  and 
the  theory  in  Billingsley  directly  applies.  In  this  paper,  the  theory  is  ap¬ 
plied  to  make  inferences  about  how  arrival  and  service  rates  vary  with  the 
number  in  the  system. 

The  birth  and  death  processes  considered  here  are  restricted  so  as  to  be 
ergodic  and  metrically  transitive.  A  likelihood  function  is  developed  based  on 
observing  a  system  continuously  over  some  interval  of  length  T  .  Relations 
which  are  needed  to  work  out  results  for  particular  models  are  developed,  and 
limit  theorems  from  Billingsley  are  presented.  Estimators  and  hypothesis  tests 
are  derived  for  three  simple  models,  as  well  as  their  asymptotic  properties. 

All  the  models  considered  will  have  only  a  finite  number  of  parameters. 

Of  particular  interest  are  composite  hypotheses  such  as  that  the  service 
rate  does  not  vary  with  the  number  in  the  queue.  A  conrposite  hypothesis  can 
be  tested  by  first  nesting  it  in  a  more  general  hypothesis.  Methods  for  doing 
this  are  discussed  at  several  points  in  the  paper.  A  numerical  example 


illustrating  these  results  is  then  presented,  which  includes  the  testing  of  both 
true  and  false  composite  hypotheses. 

More  general  models  for  the  variation  of  arrival  and  service  rates  with 
the  number  in  the  system  are  considered  in  an  Appendix,  as  well  as  other  ideas 
for  further  research. 


II.  Definition  and  Characteristics  of  the  Process 

A  queueing  process  with  Poisson  arrivals  and  exponential  service  is  a 

birth  and  death  Process.  If  the  number  of  customers  waiting  in  the  queue  plus 

the  number  in  service  is  j  ,  we  say  the  system  or  process  is  in  state  E,  . 

J 

When  the  system  is  in  E  arrival  and  service  rates  will  be  denoted  by  X  and 

"  J 

(X.  respectively.  Hence,  X  and  |j.  may  be  functions  of  E  but  given  E  , 

**  J  j 

they  do  not  vary  with  time.  (They  are  time  homogeneous.  )  If  a  customer  arrives 


while  the  system  is  in  state  E.  ,  we  have  transition  E. 

J  J 


Ej+1  .  Similarly, 


if  a  service  is  completed  while  in  state  E.  ,  we  have  transition  E  — «■  E.  . 

J  J  J  "1 

In  all  models,  the  variation  of  Xj  and  with  Ej  will  be  expressed  in 

terms  of  a  finite  number  of  parameters,  6'  =  (©^©g,...,©  )  . *  We  will  assume 
the  following  throughout  this  paper; 


1. 

a  0 

2. 

ti. 

>° 

for 

J  B  1,2,... 

3. 

*J 

>o 

for 

all  Integral 

J  <M  , 

where 

M  may  be 

4. 

*J 

a  0 

for 

J  =  M  ,  M  + 

1 

if  M 

is  finJ  te 

5. 

The 

mean  recurrence  time 

for  each 

state 

is  finite. 

Assumptions  1-4  prevent  the  existence  of  two  disjoint  non-null  irreducible  sets  of  state 
That  is,  the  process  must  be  metrically  transitive.  Should  metric  transitivity 


In  this  paper,  ©  denotes  a  column  vector,  while  its  transpose  ©'  ,  denotes 
a  row  vector. 
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not  hold,  the  theory  will  apply  to  each  irreducible  subset  of  states  considered 
separately.  With  the  addition  of  Assumption  5,  the  process  will  be  ergodic. 
Conditions  on  the  {'X.}  and  {|J..}  for  Assumption  5  to  hold  are  known,  but  of 
little  importance  here  since  the  dimensionality  of  the  parameter  space  is  not 
affected,  and  the  mean  recurrence  times  will  obviously  be  finite  for  most  ap¬ 
plications.  In  order  not  to  violate  Assumption  5,  we  must  choose  models  with 
the  property  that  for  some  9’  there  exists  a  neighborhood  such  that  any  O' 
belonging  to  that  neighborhood  defines  a  process  which  is  ergodic  and  metrically 
transitive.  A  queueing  model  which  is  not  ergodic  is  the  following: 

2 

Xj  =  a(j )  J  =  0,1,. . .  a,b  >  0 

!Lj  “  h j  j  a  1,2,... 

where  the  unknown  parameters  are  O'  *  (a,b)  . 

Ergodicity  and  metric  transitivity  quarantee  that  large  sample  theory  is 

meaningful .  Suppose  the  system  is  in  state  Ej  at  t  =  0  .  We  represent  the 

probability  that  the  system  is  in  state  E^  at  t  by  PfE^ft)  |  Ej(0)]  =»  Pjk(t) 

Under  the  above  assumptions,  Lim  P ,,  (t )  =  P  >  0  ,  independent  of  J  such 

-U.00  Jk  K  ^ 

that  (P,  }  is  the  unique  stationary  distribution  for  the  process,  2  p  =  1  . 

k=0  * 

Other  properties  of  birth  and  death  processes  are  discussed  in  Feller  [5]. 

III.  Observing  the  System 

It  will  be  assumed  here  and  in  what  follows  that  the  system  is  observed 
continuously  over  some  time  interval  T  .  During  T  ,  the  times  at  which 
customers  arrive  and  services  are  completed  are  observed.  This  is  equivalent 
to  considering  the  state  of  the  system  as  a  step  function  (not  monotone)  over 
time.  Suppose  one  starts  to  observe  the  system  at  t^  ~  0  and  finds  the 
system  in  state  E^  where  j  =  xi  .  Let  the  first  transition,  Ej  -* 
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(an  arrival  or  completion  of  service)  occur  at  tg  ,  where  k  -  Xg  .  Let  the 
sequence  of  transitions  after  tg  occur  at  ty . . .  ,tn+^  >  tQ+g, . .  .  and  let 
the  number  in  the  system  instaneously  after  t^  be  .  If  t  ^  £  T  and 

tn+2  >  T  »  then  one  ***  the  se<luenoe  (X2J  *2 ) » *  *  ‘ »  (Xn+1' *n+l  ^  *  where 

n  is  the  number  of  transitions  observed  during  T  . 

IV.  Constructing  the  Likelihood  Function 

In  order  to  construct  the  likelihood  function  for  the  above  sequence,  we 
will  first  consider  a  single  transition,  x^  — *  Xi+1  *  Note  that  because 
transitions  are  single  step,  |  *i+1  "  |  ■  1  •  Let  x^  »  J  and  let 

■  ti+1  -  ,  the  time  the  system  remains  in  state  Ej  ,  upon  entering  at 

t±  . 

That  portion  of  the  likelihood  function  due  to  a  single  transition  is  made 
up  of  two  components,  the  density  function  of  time  in  state  Ej  evaluated  at 
;  and  given  ,  the  probability  that  the  transition  at  ti+1  is  the  one 
observed  (an  arrival  or  completion  of  a  service).  Given  x^  ,  the  probability 
that  >  T  is  simply  the  probability  that  both  the  next  arrival  and  next 
service  completion  occur  after  t  .  Since  these  events  are  independent, 

-M-.T 

P[t±  >  t]  =  (e  )(e  )  *  exp 

Hence,  the  density  function  for  is: 

f('r±)  -  (>±  +  \x±)  exp  -(7^+p^)^ 

Note  that  we  are  using  \  and  p  to  refer  to  the  state  of  the  system  E 

1  1 

in  ft  .t  ]  ,  and  not  to  any  particular  transition.  Strictly,  one  should 

\  j r 

use  1  and  p  ,  but  this  notation  is  needlessly  cumbersome. 
xi  xi 
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It  Is  not  difficult  to  show  that  P[x1+1  =  x^l  |  T^]  »  +  n  )  , 

Independent  to  T±  .  Equivalently,  the  probability  that  the  next  transition 
is  due  to  a  service  is  . 

Hence,  we  obtain  the  portion  of  the  likelihood  function  due  to  a  single 
transition,  dF±  «*  dF4(xi+1  ,  |  x^  ,  0 )  ,  in  the  following  form: 

-  l\/(\  +  |d± )  ]  (\  +  M-±  )exp  -(>>i-Hi1)T1 

.  exp  -(>^4^)-^  for  x±+1  -  x±  +  1 

-  M-±  exp  -(>>1-hj.1)t1  for  x±+1  -  x±  -  1  (l) 

The  correct  form  depends  on  what  was  actual 3 y  observed  (an  arrival  or  a  service). 

Since  the  future  behavior  of  a  Markov  Process  depends  only  on  the  state  of 

the  system  when  last  observed,  the  portion  of  the  likelihood  function  due  to 

a  single  transition,  conditioned  on  the  state  of  the  system,  is  independent 

of  the  past  history  of  the  system.  Hence,  the  likelihood  function  of  the  entire 

sequency  of  n  transitions  is  simply  the  product  of  terms  such  as  (l).  In 

addition,  there  is  the  probability  that  the  system  started  in  state  E  .  If 

*1 

the  system  is  initially  observed  at  a  random  moment  along  the  real  line, 

P[x-L]  a  P^  ,  the  stationary  probability.  However,  for  large  sample  theory, 
the  assumption  of  random  entry  is  not  critical.  Further,  the  portion  of  the 
likelihood  function  due  to  transition  n  +  1  occurring  after  T  is 
P[Tn+1  >  T  -  tR+1]  =  exp  "(\1+i'HJ-n+1)(T-'fcn+1)  •  Hence,  the  likelihood  function 
for  the  entire  set  of  observations,  Ln (9)  ,  is  the  following: 

Ln(i)  =  PtxJ  exp  -(Vl^n+l)(T-tn+l)  “  ^i  •  (2> 

issi 

It  is  convenient  to  deal  with  the  log  of  the  likelihood  function  lnL^C©)  . 
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As  T  — ►  oo  ,  n  — ►  oo  with  probability  one,  and  for  large  sample  theory,  the 
first  term  in  lnLQ(0)  ,  lnPt^]  ,  can  be  ignored.  When  this  is  done, 
lnLn(0)  can  be  put  in  the  following -very  useful  fdrm: 

\ 

00  oo  00 

inLn(e)  .  £  ♦  I  Vt  -  I  r,(*j  *,  Hj)  ,  (3) 


J-0 


J«1 


j— 0 


where  u.  is  the  ntmber  of  transitions  E. 

J  J 


Ej+1  (up)  ,  dj  is  the  number 


of  transitions  E.  — ♦  E  (down)  ,  and  y .  is  the  total  time  spent  in  E 
J  J “l  j  J 

during  T  . 


V.  Useful  Relations  for  Later  Application 

In  what  follows  A,  and  n  will  be  allowed  to  vary  with  the  number  in 
J  J 

the  system  in  a  way  that  can  be  functionally  described  in  terms  of  a  finite 
number  of  independent  parameters.  Let  0'  -  (Q  . . . . ,  Q,)  be  the  vector 

1  K 

representation  of  the  parameters  in  Euclidean  k-space.  We  have  in  general 

Tvj  =>  hj(®)  >  M-j  =  gj  (6) ,  suitably  constrained  as  in  Section  II. 

In  addition  to  the  stationary  distribution,  [P.)  ,  we  will  need  the 

J 

distribution  C^Tj ) ,  where  7Tj  is  the  probability  that  a  randomly  selected 

transition  was  out  of  state  E.  •  With  greater  precision,  define 

7T  =  Lim  (u  +  d  )/n  ,  which,  for  processes  of  the  type  considered  here  exists 

J  T-*  CO  J  J 

with  probability  one.  For  an  ergodi:  process,  a  limit  such  as  7r.  is  inde- 

J 

pendent  of  the  starting  state,  or  equivalently,  of  any  distribution  of  starting 

states.  In  evaluating  this  limit,  however,  it  is  conveneint  to  assume 

starting  conditions  are  determined  by  the  stationary  distribution,  (Pj)  . 

When  this  is  done,  we  obtain  for  any  interval,  T  ,  E[(u  +  d  )  j  T]  =»  (X  +  u  )PT 
00  J  J  J  J  J 

and  E[n  |  T]  =  I  (X,  +  (j.  )P  T  .  It  is  now  clear  that 
j-0  0  J  0 
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* 

Under  suitable  regularity  conditions,  and 
where  Cl  denotes  the  entire  sample  space, 
are  obtained : 


from  £ dFi(x1+1,T1  |  x±,9)  -  E(l) 
the  following  well-known  results 


1 


E(GU)  -  0  (6) 

and 

E(G  )  -  -a  .  (7) 

'  uv  uv 

From  (7),  we  see  that  the  variance -covariance  matrix  of  Gu  ,  E(GG‘ )  is: 

e(gg’)  =  a(e)  -  (ouv)  =  (-e(guv))  •  (8) 

VI.  Limit  Theorems  from  Billingsley  [l] 

THEOREM  1:  [1]  The  asymptotic  distribution  of  n _1/2  as  T  — *  00  Is  the 
multivariate  normal  with  mean  9  and  variance -covariance  matrix  o(9)  .  We 
write : 

n"l/2  §a  ^  N(9  ,  0(9))  . 

A  A  /V  . 

Next  consider  the  maximum  likelihood  estimator,  9’  a  (9  ,  .....9  .  )  ,  found 

— n  m.  nk 

by  setting  ■  0  .  Assume  that  for  large  n  the  estimator  is  unique. 

THEOREM  2:  [1]  The  maximum  likelihood  estimator,  9  ,  converges  in  probability 
as  T  — »  00  to  the  true  vector  of  the  parameters,  9  ,  such  that : 

n1//£(9  -  9 )  N(0  ,  0_1(e)) 


Specifically,  conditions  on  dFj.  which  permit  the  interchange  of  the  order 
of  the  operations  of  differentiation  and  integration.  See  Cramer  [4,  p.  67 
and  501],  or  Billingsley  [1,  p.  5  and  4l]. 
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o 

If  we  hypothesize  that  ©  is  some  particular  vector,  ©  say,  such  that 
the  probability  distributions  of  all  random  variables  in  the  process  are 
completely  specified,  then  we  have  a  simple  hypothesis.  Given  :  ©  ■  ©°  , 
the  vector  ©  is  restricted  to  a  single  point,  of  dimensionality  zero,  in 
the  k-dimensional  parameter  space.  A  general  way  to  test  such  a  hypothesis  is 
with  the  likelihood  ratio. 


or  equivalently,  with  minus  twice  the  log  of  this  ratio: 

2 [Max  lnLn(©)  -  lnLn(©°)]  .  (9) 

If  the  value  obtained  from  (9)  exceeds  some  critical  value,  we  reject  Hq  .  We 
would  like  to  know  the  asymptotic  distribution  of  (9)  when  HQ  is  true,  and 
also  under  alternative  hypotheses. 

THEOREM  3:  [1]  If  HQ  :  ©  ■  ©°  is  true,  then  minus  twice  the  log  of  the 
likelihood  ratio  is  asymptotically  distributed  as  chi  square  with  k  degrees 
of  freedom: 

2 [Max  lnL^©)  -  lnLn(©°)] 

Suppose  that  Hq  is  false,  and  some  alternative,  :  ©  ■  61  ,  is  true. 

Define  the  vector  h  as  follows : 

9°  -  Q1  +  n"1//2  h  or  h  =  n1/^©0  -  01)  . 

on 

How  define  the  quadratic  form,  A  =  h'o(©  )h  ,  where  A  is  the  non-centrality 
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parameter  in  the  non-central  X  distribution  below. 

THEOREM  4:  [1]  If  Hq  :  9  ■  6°  is  false,  but  :  ©  ■  91  is  true,  then; 

2 [Max  lnLn(S)  -  lnLn(£°)J  X^A  . 

Suppose  one  imposes  restrictions  on  the  components  of  9  ,  3uch  that  some 
of  them  are  equal,  or  are  known,  or  are  related  to  one  another  in  such  a  way 
that  given  an  appropriate  proper  subset  of  the  {G^}  ,  the  rest  can  be  solved 
for  uniquely.  When  such  is  the  case,  one  can  represent  the  system  in  terms  of 
a  vector  of  k '  independent  parameters,  £  ,  such  that  k '  <  k  .  A  model 
incorporating  these  new  restrictions  is  said  to  be  nested  in  the  original 
model,  and  is  represented  by  0c0 .  The  components  of  £  can  be  estimated 
and  tested  by  methods  used  for  9  above.  In  addition,  the  hypothesis  that 
the  restricted  model  is  true  under  the  assumption  that  the  original  model  is 
true  can  also  be  tested. 

THEOREM  5:  [1]  If  the  restricted  model  is  true,  and  the  hypothesis  that 
^  a  gt°  is  also  true,  then: 


2 [Max  lnL  (9°)  -  lnL  (6°)]  =*  2 [Max  InL  (9)  -Max  InL  (£)] 
o  n_  n_  a  n  —  d  n 


9  “  £ 

+  2  [Max  lnLn(£)  -  InLjg0)] 


O  JO 

where  G  is  determined  by  J  ,  and 


2  [Max  lnLn(9)-  lnLn(0°)]  ^  xjj 
2[Max  lnLn(£)-  lnLn(g°)]  X^, 

l 

2[Max  lnL  (9)  -  Max  lnL  (£)]  X^_k, 

9  £ 
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where  the  last  two  expressions  are  asymptotically  independent. 

In  the  remainder  of  this  paper,  the  theory  is  applied  to  birth  and  death 
problems  of  interest. 

VII.  Application  of  the  Theory  of  Particular  Models 

MODEL  1:  Xj  >  0  ,  J  -  0,1, . . .  ,M-1  XM  -  ° 

(Xj  >0  ,  J  *  1)2,  •  •  •  ,M  |1q  **  0 

As  mentioned  earlier,  restrictions  must  be  placed  on  models  we  consider 
to  explain  the  variation  of  the  (Xj)  and  the  {(Xj)  With  J  so  that  there 
are  only  a  finite  number  of  unknown  parameters.  Model  1  accomplishes  this  by 
truncation,  so  that  the  number  of  the  system  cannot  exceed  M  .  There  are 
2M  unknown  parameters  (Xj)  and  (pLj)  in  all,  which  are  required  only  to  be 
finite  and  positive. 

The  maximum  likelihood  estimators  are  determined  from  Equation  (3),  which 
for  this  model  has  the  form: 

M-l  M  M 

inLn(-) "  I  ujlnV  I  djlnp'j  I  rj(V^}  (10) 

j=0  J-l  j=o 


Differentiating  with  respect  to  each  of  the  unknown  parameters,  we  get: 


3LnL 

n 


and  therefore 


for  j  =>  0,1,..., M-l 


(11) 
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31  nL  d. 

J.  .  r  -  0 

d(1J  J 

and  therefore 


for  j  m  1,2, ...  ,H  .  (12) 


There  is  nothing  particularly  surprising  about  (11)  and  (12).  The  estimate  of 

7V  is  the  number  of  arrivals  while  in  state  E  divided  by  the  total  time 
<3  J 

spent  in  state  E.  .  The  analysis  does  not  stop  here,  however. 

J 

For  the  variance -covariance  matrix  of  the  estimators,  we  begin  with 
Equation  (l)  and  Gu  -  3/d«u  lndF±  ,  and  given  -  j  obtain: 


-  T. 


for  x1+1  ■  j  +  1 


for  •  J  -  1 


and  all  other  *  0  . 

The  procedure  to  determine  the  variance-covariance  matrix  of  the  (pu)  > 
o(Q)  ,  will  be  first  to  determine  0(9  |  for  (7).  Then,  from 

E^  |  x^)  =  E^)  -  0  for  all  u  ,  E(o(e  |  x± ) )  =  o(6)  .  We  continue  to 
condition  on  x^  •  j  : 


Vj 


0 


f°r  Xi+1  =  J  +  1 


0 


for  x,  .  *  j  -  1 
i+1 


and  all  other  Guy  =  0  . 
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As  shown  earlier 


,  P(xi+i  -  J  +  1  |  -  J)  -  +  M-j)  • 


Now 


V(G. 


^ 1  *»  ■ 11  ”  -E(v,  ■ J)  ■  (1A*>  *r*7  ■ 


for  J  m  0,...,M-1 


V(G4i  I  x±  -  J)  -  -E(G^  I  X±  «.  J)  -  (1/n*)  -  HjftjVuTJ 


J+M-j  HjVAj+Uj- 


for  j  -  1, , . .  ,M  ,  and  Cov(Gu,Gy  |  x±  -  J)  -  -E(Guv  |  x±  -  J)  -  0  for  all  other 

(u,v)  pairs.  Also,  V(G^  I  ^  ft  J)  *  V(G^  I  ^  J)  •  0  . 

J  J 

We  now  obtain: 

M  n 

T<V  .  E(V(0^  I  -  *»  -  X  V(0X,  I  *1  ■  kK  ■  iv,) 


k»0 


for  J  -0,1,...,M-1 


7t 


V(G  )  -  E(V(G  |  x,  =  k)}  =  — rr — j.  -,;— y 

Vj  V)  1  +  *V 


for  J  =  1,2,...,M  ,  and  Cov(Gu,Gy)  -  0  .  Hence,  using  Equation  (4), 
P, 


v(V ' 


for  J  =  0,1, . . . ,M-1 


v(V  -  5^ 


for  3  m  1,2, ... ,M 


Since  a.1 3  the  covariances  are  zero,  o(0)  is  a  diagonal  matrix.  Let 
G'  -  (G^,...,G^  ,0^,...,°^)  ,  then  0(©)  -  (l/2B) (^j)  where 
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9 


l'  -  (PqAq,  . .  *  • '  *  '  W  •  NoW  -  (°UV>  -  2R(&i/ii) 

and  letting  ©'  -  (XQ,  ••  •»V-1,^V*  •  *  >  we  by  Theorem  2: 

nl/2(©  ‘  ®)  *  N(0  ,  (T^©))  . 

Since  c-1(©)  is  a  diagonal  matrix,  the  maximum  likelihood  estimators  are 

asymptotically  independent.  If  we  let  V  (X)  be  the  asymptotic  variance 
“  as 

(the  variance  of  the  asymptotic  distribution)  of  X  ,  then: 

2^  R 

V»'1/26.1  -  V5  "  “p^*  f°r  J  -  0,1,..., M-l  (U) 

j 

^  R 

V«sln’1^J  ■  ■  ~T~  f°r  i  ■  1,2,..., M  (14) 

U 

Suppose  we  hypothesize  that  9  is  some  particular  vector, 

©°  -  (?£,..., x£_1,m£,...,|oJJ)'  •  If  HQ  :  ©  -  ©°  is  true,  then  by  Theorem  3, 

2[Max  lnLn(9)  -  lnLn(©°)]  , 

9 

where  2M  is  the  number  of  parameters  (components)  of  0  .  For  this  model, 
we  obtain: 

M-i  A 

2[Max  lnLj©)  -  lnLn(0°)]  -  +  L  } 

2-  J-0  J-l 

+  I  ri^°i  +  ^V  '  n 

j-0 

which  is  computable  from  the  sample.  We  determine  a  critical  region  such  that 
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2  o 

P(x^  c)  =  a  ,  then  reject  Hq  :  0  ■  9  if  the  number  obtained  from  (15) 

exceeds  c  . 

If  Hq  is  false  such  that  ©  =  91  say,  then  the  power  of  the  test  can 

be  determined  as  a  function  of  9°  and  01  .  For  h  -  n1/2(9°  -  01)  ,  the 
2  1 

quantity  A  =  h'a(©  )h  is  easily  computed: 


a2  11 

A  “  2E 


.  ^ 


PM-1^M-1  '  *M-l)  P1 W  '  M^) 

+.  ..+  - J -  +  - J - h...  + 

*m-i  ^ 


T 


% 


(16) 


From  Theorem  ^  we  have : 


la^(£)  -  InL^g0)]  £. 

Given  0°  and  01  ,  A  can  be  computed  from  (16),  and  then  the  power  of  the 
test  determined  as  ^  >  c  )  . 

Nested  within  this  model  is  the  simple  queue  limit  model  where  it  is  assumed 
that  Xj  and  (ij  do  not  vary  with  the  number  in  the  system.  Suppose  one  owns 
a  gas  station  and  has  an  opportunity  to  enlarge  his  facility  by  buying  adjacent 
land.  How  much  additional  business  will  he  get?  Does  he  lose  business  only 
when  the  facility  i^  full  (state  E^)  ,  or  are  customers  also  discouraged  by 
queue  sizes  less  than  M  -  1  ?  Without  assuming  any  relationship  among  the  u,j  , 
the  hypothesis  that  =>  =  . . .  =  ^  can  tested,  where  M  ife  the 

maximum  number  of  cars  the  facility  can  hold.  Let  <£'  =  . . .  ,10^)  . 

From  (10),  the  likelihood  of  the  restricted  model  is: 

M-l  M 

djlDM'j  "  +  ^j1 

J*»0  Jsl  j=0 


in^(^)  "  ln*  X  + 
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and 


M-l 

^  lnLn(^)  -  ^  Uj  A  -  (T  -  rM)  -  0 

J-0 

therefore 


M-l 


which  is  easily  computable  from  the  sample. 

MODEL  2:  ?Vj  -  ah(j)  J  -  0,1,...  and  g(0)  »  0 

M-j  “  ^g(j) 
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where  a  and  b  are  unknown  parameters,  and  h(  j  )  and  g(j)  are  assumed 
to  be  arbitrary  but  known  functions  of  ,3  . 

The  advantage  of  the  model  is  that  it  is  not  necessary  to  impose  a  queue 
limit  in  order  to  keep  the  number  of  parameters  finite.  Furthermore,  informa¬ 
tion  about  a  and  b  is  obtained  throughout  T  ,  rather  than  only  when  the 
system  is  in  certain  states.  Hence,  statistical  inferences  will  be  more 
precise.  However,  in  order  to  use  such  a  model,  considerable  information 
about  the  behavior  of  the  system  mu3t  already  be  at  hand.  When  a  queue  limit 
is  imposed  by  the  nature  of  the  system,  then  Model  1  is  more  general  than 
Model  2. 

For  Model  2,  (3)  takes  on  the  following  form; 

00  08  00 

lnLn(S)  -  ^  Ujln  ah( j )  + 

J-0  J«1  <3*0 

from  which  we  obtain: 


^  djln  bg(J )  -  ^  Tj [ah( J )  +  bg(j)] 


(18) 


00  /  00 

Note  that  Lim  Z  u.  /  Z  d.  *  1  with  probability  one,  for  in  ergodic  queues, 
J-0  V  J=1  0 

whatever  goes  up  must  come  down.'  For  large  samples,  we  can  therefore  replace 

00  oo 

Z  u.  and  Z  d.  by  n/2  . 

J-0  J  J-l  0 

We  now  obtain  the  variance -covariance  matrix,  starting  again  with  the 
likelihood  of  a  single  transition: 
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dF(xi+l,Ti  I  xi  “  ^»£)  “  exP  -(ah(j)+bg(j))T±  if  x±+1  -  j  +  1 

-  t>g(J)  exp  -(ahCjJ+bgO))^ 


if  xi+l  -  i  -  1 


and  given  »  j  , 


Ga 

u 

9)  |H 

-  hO)^ 

,  Gb  -  -  e(2)\ 

if 

xi+l 

■  J 

+  1 

G 

a 

- 

-  b(J)T± 

>  Gb  "  b  ■  g(J)Ti 

if 

Xi+1 

-  J 

-  1 

Gaa 

-  - 

1 

T 

a 

Gbb  - 0 

if 

X.  , 

i+1 

-  J 

+  1 

G 

.  0 

G..  =  - 

if 

X.  . 

=  J 

-  1 

aa 

bb  b2 

i+1 

°ab  ■  0  ' 


Hence : 

V(Ga  I  x±  -  j)  -  (1/a2)  bgO)] 

v(Gb  I  xi  "  J)  "  O-A  )  ah(j)S+Jlg(j) 

Cov  (Ga,Gb  |  x.±  m  j)  .0  and 


bg(j)J 


V(Go 


oo  w 

a>-  I  ,<0JIr)>V  I 


h(J>n- 


J-0 


J-0 


V(Gb)-  £  V(Gb  |  x±  -  Jkj  -  £ 


a[ah(j)  +  bg ( J ) ] 
g(J>7T1 

b[ah(j)  +  bg(j)] 


j-0 


J-0 


How  from  (4),  7T  -  P,[ah(j)  +  bg(j)]/2R  ,  and  we .  obtain  the  simple  result: 
J  u 
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1 


v(0a) 


00 


h(J)Pj  - 


Cov  (Ga,Gb)  -  0  . 

For  0'  a  (a,b)  , 


again  a  diagonal  matrix. 

Thus  for  n1//2(0  -  9 ) '  -  n1//2(a  -  a  ,  b  -  b)  ,  n1/2^  -  0)  N(0  ,  c_1(0)) 

where 


Note  that  the  form  of  this  result  is  independent  of  h(,j)  and  g(j)  .  However, 
this  does  not  imply  that  the  variances  are  independent  of  h(j)  and  g(j),  since 
for  given  T  ,  the  distribution  of  n  depends  on  them. 

Hypotheses  about  a  and  b  can  also  be  tested  by  the  same  methods  used 
in  Model  1.  If  0°  ■  (a°,b°)'  is  hypothesized,  but  01  a  (a^b1) '  is  true, 
such  that  h  a  '/n(0°  -  O^)  ,  then: 
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2  [Max  lnLn  (6)  -  lnLn(0O)]  —  x|>A 


Special  Cases  of  Model  2  Which  are  of  Independent  Interest;  Included  in  Model  2 
are  many  common  queueing  models.  Some  of  these  Eire  listed  below: 

Model  2a:  The  infinite  channel  model,  Xj  =  X  ,  Pj  =  JM-  >  for  «5  =  0,1,...  . 

This  is  simply  Model  2  with  X  =  a,h(j)  =  1  ,  b  =  p  ,  g(j)  =  j  . 

Model  2b;  The  C-channel  model,  Xj  “  X  >  M-j  ”  JM-  Tor  3  *  0,1,...,C-1 
Pj  =  qjL  for  j  =  C,C+1, ...  .  This  also  is  Model  2  with  X  »  a,h(j)  -  1  , 
b  =  p  ,  and  g(j)  a  j  for  j<C-l,g(j)=C  for  j  ^  C  . 

Model  2c;  The  machine-repair  model,  single  repairman  case.  Suppose  there  are 
M  machines,  each  of  which  break  down  at  rate  X  when  operating.  The  machines 
are  in  operation  continuously  except  when  they  break  down.  The  repairman  works 
on  one  machine  at  a  time,  at  rate  p  .  Let  the  state  of  the  system,  Ej  ,  be 
defined  as  the  number  of  working  machines.  Then  we  have:  Xj  =  jX  for  j  {  M  , 
p_  8  p  for  j  £  M-l  ,  -  0  .  Note  that  a  service  corresponds  to  a  transi¬ 
tion  Ej  -*  Ej+i_  while  an  arrival  results  in  a  transition  Ej  -*  •  Thus 

we  have  a  truncated  C-channel  model  with  the  roles  of  X j  M-j  interchanged. 

Note  also  that  this  model  is  nested  in  Model  1  . 

Model  2d;  Model  2c  with  N  reparimen  such  that  N  <  M  ,  and  each  repairman 

works  alone  at  rate  p  .  For  this  case:  Xj  *  jX  Tor  0  <  j  <  M  ,  Pj  =  (M-l)p 

for  M-N<j^M,  p  =  Np  for  O  <  j  ^  M  -  N  .  This  model  is  also  nested 

in  Model  1,  and  represents  the  hypothesis  that  the  rate  repairmen  work  is  not 
affected  by  the  number  of  jobs  waiting.  Assuming  Model  1  correctly  represents 
the  system,  this  hypothesis  can  be  tested  by  the  methods  presented  earlier. 
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Rejection  of  the  hypothesis  leads  to  the  conclusion  that  the  number  of  Jobs 

waiting  does  affect  the  performance  of  repairmen.  It  would  be  more 

efficient,  however,  to  nest  this  model  in  a  model  more  restrictive  than  Model  1, 

assuming  that  the  new  model  is  true.  One  such  model  is:  A.  *»  JA  ,  and 

J 

4j  -  0  1  J  1  M  . 

Comparing  Model  1  and  Model  2,  we  see  that  Model  1  is  more  general,  but  if 
the  queue  limit  must  be  made  large  to  be  realistic,  the  number  of  parameters 
becomes  correspondingly  large.  Under  these  conditions,  tests  lose  power,  and 
particularly  for  large  J  ,  the  variances  of  estimators  are  apt  to  be  large.  ' 
On  the  other  hand.  Model  2  is  rather  restrictive  and  requires  a  lot  of  a  priori 
information  to  be  useful.  A  compromise  is  suggested.  Model  3  is  such  a 
compromise. 

MODEL  3:  Aj  >  0  for  J  =  0,...,M-1  ,  p.j  >  0  for  J  “  1,...,M-1 

Aj  =  *  ^J  "  f°r  J  ^  M 

For  situations  where  large  queues  are  possible,  but  the  service  facility 

performs  well,  relative  to  the  load,  such  that  the  average  queue  length  is  small, 

this  model  will  be  quite  useful.  The  A.  and  p.  can  be  individually  esti- 

J  J 

mated  for  small  J  ,  where  the  system  is  most  of  the  time,  and  lumped  together 
for  large  j  ,  without  serious  consequences. 

From  (3),  we  obtain  for  this  model: 

Aj  =  u./Yj ,  =  yTj  (25) 

as  before,  and 


(26) 
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(27) 


oo  /  oo 

s- 1  VI  rJg<j) 

j-M  /  J-M 

In  obtaining  the  variance-covariance  matrix,  all  the  covariances  will 
again  be  zero.  The  only  part  that  is  at  all  new  is  the  determination  of 
V(Qa)  and  VfG^)  .  If  x^  »  J  <  M  ,  G  and  are  identically  zero.  Hence 

we  need  only  consider  -  J  ^  M  . 

Given  x±  -  J  £  M 

Ga  “  a  '  h(J>Ti  »  Gb  “  _s(J)Ti  if  xi+i  -  J  +  1 

Ga  -  -hG))^  Gb  “  b  ‘  g(J)Ti  if  xi+i  “  J  ‘  1  • 

It  is  now  clear  that : 
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Many  intereting  possibilities  are  uncovered  by  noting  that  Model  2  is 
nested  in  Model  3.  Thus,  one  can  test  the  hypothesis  that  a  particular  h(j) 
and  or  g(j )  applied  to  the  lower  states  only.  The  number  of  states  included 
in  such  a  test  is  determined  by  M  .  Further,  there  is  no  need  to  require 
that  M  be  the  same  for  both  the  C>v j }  and  the  (p,j)  • 


VIII.  An  Example  Using  Data  Generated  by  Simulation 

We  would  like  to  know  how  well  the  -theory,  when  applied  to  models  such  as 
those  in  section  VII,  can  discriminate  between  those  which  reflect  the  structure 
of  a  system  and  those  which  do  not.  An  analysis  of  data  generated  by  simulating 
a  model  of  known  structure  will  be  of  value  here. 

For  the  simulation,  the  following  model  was  used: 


-  4/(J  +  1) 

M-o  -  0 

Ul  -1 

“  2 


j  **  0, 1, , » * 


2,3,, 


Results  of  the  simulation  (along  with  F  (true)  for  comparison)  were  as  follows: 


j 

rJ 

P^true) 

0 

44 

0 

10.54 

.065 

.073 

1 

93 

44 

50.76 

.312 

.291 

2 

63 

93 

39.33 

.24 2 

.291 

3 

31 

63 

32.34 

.198 

.193 

4 

13 

31 

15.90 

.098 

.097 

5 

5 

13 

10.66 

.0  66 

.039 

6 

1 

5 

2.54 

.016 

.013 

7 

_ 0 

1 

_ sli 

.004 

.004 

z 

250 

250 

162.80 

1.001 

1.001 
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Note  that  u  »  d.  .  This  relationship  must  hold  for  all  states 
J  J+X 

except  initial  and  final,  and,  by  chance,  did  for  them  too  in  this  simulation. 

00 

It  is  not  difficult  to  show  that  2  P.  <  .02  P.  and  hence  is 

J-8  J  0 

negligible.  We  also  compute  the  true  throughput,  R  -  +  p.2(l  -  PQ  -  P.^)  - 

-1.56  . 

Assume  for  concreteness  that  we  have  a  2-channel  facility  for  which  it  has 

been  observed  that  rarely  are  there  more  than  4  customers  in  the  system  at  any 

time.  (it  might  be  a  gas  station,  a  bank,  or  even  a  bordello.)  We  would  like 

to  know  how  varies  with  J  .  Specifically,  for  small  j  ,  does  it  decline 

J 

with  j  or  remain  constant?  The  economics  of  enlarging  the  facility  are 
closely  tied  to  this  relationship. 

The  model  we  will  assume  in  the  analysis  is  the  following: 


Model  :  Hq  :  ^  ® 

7\j  «•  a/(j  +  1) 

M.0  =■  0 

^  =  (1 


J  =  5,6,, 


J  =  > 


Note  that  there  are  7  unknown  parameters. 

Estimation  of  parameters  and  their  coefficients  of  variation 


true  9 

est.  eq. 

A 

9 

Vas  eq. 

Vas(®  "  9) 

a/[i 

4.00 

(25) 

4.17 

(15) 

.342 

.146 

2.00 

(25) 

1.83 

(15) 

.0428 

.103 

1.55 

(25) 

1.60 

(15) 

.0285 

.126 

1.00 

(25) 

0.96 

(15) 

.0321 

.179 

0.80 

(25) 

0.82 

(15) 

.0499 

.280 

4.00 

(26) 

2.68 

(28) 

2.81 

.418 

1.00 

(20) 

0.99 

(21) 

.0040 

.064 
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Let  us  denote  the  estimate  of  a  in  this  table  by  ■  2.68  . 

Next  consider  another  model. 

Model  Mg  :  HQ  :  -  a/(j  +1)  j  -  0,1,... 

^0  "  0 
M-!  "  M- 

M-j  "  ^  J  "  S»5>*“ 

(this  assumption  is  correct)  and  we  obtain: 

-  4.00 

Vas^2  "  a>  “  *064 

(oAOfc  =*  .064  -  (a/|J.)^  as  well. 

2  J 

We  will  now  test  whether  model  Mg  is  correct  noting  that  it  is  nested  in  Mj^  . 
There  were  7  unknown  parameters  under  the  original  assumption,  while  there  are 
two,  •  (a,n)  under  hypothesis  Mg  .  Hence,  from  Theorem  5  we  obtain: 

2  [Max 
6 

From  (5)  and  the  naxl.mum  likelihood  estimators,  the  L.H. S.  becomes: 


lnLn(e)  -  Max  lnLn(jZ[)]  x|_2 


if 


*0 


is  true. 


x  +  i  uJin(V(j + 1))-  x  uj  - 1  *irj/(j + 1) 


j-0 


J=5 


j=0  j=5 


-  ^  Ujln(ag/  ( j  +  D)  +  ^  agTj/(j  +  1) 

j=0  J=o 


5.52  <  Q5  =  H.OT  . 
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Hence,  we  would  accept  this  hypothesis  at  any  reasonable  a  level. 


On  the  other  hand,  let  us  test  an  alternative  hypothesis: 


Model  Mj  :  : 


unrelated  to  a  ,  with  the  same  assumptions  about 
the  {pLj}  .  (This  model  Is  not  correct.)  This  hypothesis  Is  also  nested  in 
the  original  one,  but  with  three  unknown  parameters,  -  (\a,M-)  •  We  have 


2 [Max  lnLn(©)  -  Max  lnLn(g)]  . 

A  4  4 

Using  >  «  2  u,/  2  y.  -  1.640  ,  the  L.H.S.  becomes: 

j-0  J  J-0  J 


52  >  Vos'9-1*9 

>  \,.0005  "  20 


Hence,  we  would  reject  at  any  reasonable  a  level. 

Although  it  is  obvious  from  the  raw  data  that  is  untenable,  it  should 

be  clear  that  these  methods  can  be  applied  to  more  complicated  cases  where 
inferences  by  inspection  are  hard  to  substantiate. 
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APPENDIX 

IDEAS  FOR  FURTHER  RESEARCH 


Models  of  Greater  Generality 

These  Ideas  can  be  extended  by  partitioning  the  (Ej)  Into  subsets, 
U;L,U2,...,Um  ,  such  that  if  Ej  €  u±  ,  -  a^O)  .  Similarly,  partition 

{Ej}  into  D1,Dg,...,Dn  such  that  if  E j  e  D±  ,  |i.j  «  •  If  all  the 

parameters  are  independent,  one  obtains: 


(A-l) 


(A-S) 


(A-3) 


(A-4) 


All  covariances  are  zero. 

In  all  the  models  considered  thus  far,  all  the  covariances  have  been  zero, 
so  that  o(©)  has  been  a  diagonal  matrix.  Now  we  will  investigate  why  this 
has  been  so. 

Sufficient  Conditions  for  o(9)  to  be  a  Diagonal  Matrix 

The  matrix  o(0)  is  diagonal  <=>  E(Ouv)  =  0  for  every  (u,v)  such  that 
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u  /  the  covariances  of  the  {Gu)  are  all  zero.  The  result  is  desirable  in 

that  the  variance -covariance  matrix  of  the  estimators,  c_1(©)  >  is  simple  to 
find  and  is  also  diagonal.  The  estimators,  being  already  asymptotically  normal, 
are  then  also  asymptotically  independent.  Results  can  be  stated  in  a  simple 
form. 

Now  E(Guv)  -  E[E(Guv  |  x±  =  J)}  .  Hence,  E(Guv)  -  0  if  E(&uv  |  x±  =  J)  .  0 
for  all  j  .  This  will  occur  if  dF^  |  x^  is  separable  in  the  parameters,  or 
equivalently  if  lndFi  |  x±  is  additive  in  the  parameters.  We  say  dF  |  x±  is 

separable  in  { 9^}  if  dFi  |  x^  can  be  expressed  in  the  following  form: 

^i  I  xi  "  J  ^i+l  ’  Ti  I  x±  >  eu> 

so  that  Sv(xi+1  i  Ti  I  x^  ,  ©y)  does  not  involve  u  unless  u  ■  v  .  If  this 

is  the  case,  then: 

lndFi  I  xi  -  I  lrVxi+l  '  Ti  I  xi  »  V 

u 

and  we  say  lndj^  |  x1  is  additive  in  the  parameters.  When  dF±  |  x±  is 
separable  in  the  9^  ,  then  for 

d2lndF 

u  *  v  >  Guv  I  Xi  "  0  =>E<Guv  I  xi>  "  °^E(Guv)  "  0  • 

For  all  the  models  considedered  thus  far,  arrival  and  service  parameters 
occur  together  only  in  the  exponent,  and  there  as  a  sum.  Hence,  they  are 
separable.  Further,  given  the  state  of  the  system,  at  most  one  arrival  and 
one  service  parameter  appears  in  dF^^  . 
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Not  only  have  the  likelihood  functions  been  separable,  but  the  corresponding 
models  have  been  linear.  That  is,  the  arrival  and  service  rates  have  been  linear 
functions  of  the  unknown  parameters.  The  analysis  of  variance  and  regression 
are  based  on  linear  models,  which  are  surprisingly  adaptable  in  their  ability 
to  cover  many  situations.  Further,  deviating  from  linearity  here  results  in 
considerable  algebraic  complexity.  Thus,  only  linear  models  have  been  consid¬ 
ered.  Deviating  from  separability  also  increases  complexity,  but  to  a  lesser 
degree.  The  nonseparable  linear  model  is  considered  briefly  below. 

General  Linear  Model 


Zal 


where  8Z(0)  =0  for  all  z  .  From  (2.3)  we  obtain: 


00 

V  1  | 

Y 

r — > 

00 

> 

lnLn(i)  -  L  ujln 

L  WtJ) 

+  L  djln 

L  Vz(J) 

j=0 

y=i 

j=i 

Z=1 

00 

Y 

r — > 

A 

■  I  rj 

1  ayhy(J)  + 

L  b^(J) 

j=0 

Ly=1 

z=l 

and 


(A- 5) 


SLnLjQ) 


W 

l 

j-0 


w(j) 

Y 

Z  ah  (j) 
y=i  y 


^  TTjhy(J)  «  0  for  y  -  1,2,  ...,Y 

j=o 
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&nL  (8 ) 


-n(0)  V  diSz(J)  V 
“ - -  2_,  ~Z - - - Z_,  “  0  for  z-l,2,...,Z  . 


“Sb 


j-0  Z  b  g  (j)  j-0 

z-1  z  z 


Although  the  above  equations  will  have  to  he  solved  numerically  to  obtain 
the  likelihood  estamators,  the  variance -covariance  matrix  of  the 

asymptotic  distribution- is  relatively  easy  to  find.  Below  we  find  a  (0) 

Given  =  J  , 


Ga  “  Y 
y  Z  a  h  (j) 
y-i  y  y 


MJ) 

— - VJ)T1 


if  xi+1  -  j  +  1 


hy(J)Ti 


if  xi+1  -  j  -  1 


Gb  --bz(J)T^ 


if  xi+1  =  j  +  1 


b2(J) 

~z~^ - VJ)Ti 

z  bzgz(j) 

Z=rl 


if  xi+1  =  J  -  1 


therefore 


V(Ga  1  \  =  J) 


^  ayVJ>)  l*  ayhy(J)  +  A  V.<J» 


2=1 


and  from: 


7T. 


pjU  vv(j>  *  ^(>) 


2R 
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ve  obtain 


v(g 


V 


oo 

>-4  I 


J-o 


~Y 

Z 

y-i 


hy<-»F.1 


ayhy(J) 


Similarly, 


Cov(G 


(  x  -  j) 


Y 

Z 

y«i 


yj)yj) 


Vy(J> 


Y 

Z 

yWL 


ayVJ) 


Z 
Z  b 
z=l 


zgz<J> 


therefore 


C°v(G  ,  0  ) 

y  z 


yjih.oip, 

*  »yM)) 


Similar  expressions  may  be  obtained  for  the  G,  .  Note  that  Cov(G&  ,  )  =6 

xc  y  z 

Queueing  Processes  of  Greater  Generality 

By  restricting  the  class  of  queueing  processes  in  this  paper,  it  has  been 
possible  to  apply  work  done  in  Markov  processes  to  problems  of  interest  in 
queueing.  The  author  has  also  obtained  results  for  certain  more  general 
queueing  processes,  those  which  possess  embedded  Markov  chains  [7].  One 
inport  ant  process  of  this  type  is  single  channel  with  Pols  son  arrivals  and 
arbitrary  service.  This  work  will  be  published  at  a  later  date. 
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